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ABSTRACT 

We present a new method to determine the age spread of resolved stellar populations in a starburst 
cluster. The method relies on a two-step process. In the first step, kinematic members of the cluster 
are identified based on multi-epoch astrometric monitoring. In the second step, a Bayesian analysis is 
carried out, comparing the observed photometric sequence of cluster members with sets of theoretical 
isochrones. When applying this methodology to optical and near-infrared high angular resolution 
Hubble Space Telescope (HST) and adaptive optics observations of the ^5Myr old starburst cluster 
Westerlundl and ~2Myr old starburst cluster NGC3603YC, we derive upper limits for the age 
spreads of 0.4 and 0.1 Myr, respectively. The results strongly suggest that star formation in these 
starburst clusters happened almost instantaneously. 

Subject headings: Hertzsprung-Russell and C-M diagrams — open clusters and associations: individual 
(Westerlundl, NGC3603YC) — stars: evolution — stars: formation 



1. INTRODUCTION 

Our understanding of star forma tion has progressed 
considerably in recent years (e.g., iKuiper et al.l l201ll : 
iCommerpon et al]|201lD . though the exact sequence of 
the formation of individual stars during a star formation 
event is not yet entirely understood. In particular, it is 
still unknown if low- mass stars ten d to form later than 
high-mass stars fe. g. . lKlessenll200H ) . or if high-mass stars 
form last, resulting in a rapid term ination of star forma- 
tion Ce.g.. lZinnecker fc Yorkil2007D . 

The age spread of a cluster's stellar population is a 
good indicator of the overall duration of the star for- 
mation process in the cluster. Studies of star-forming 
regions have reporte d different re sults: from a single 
age, as in NGC 410 3 (lForbeslll996D to a ge spreads of 2- 
4 Myr, e.g., i n LH95 (iDa Rio et al.ll20"Tol) . the Orion Neb- 
ula Cluster (jReggiani et al.l 1201 ID . W3 Main (|Bik et al.l 
Holl), and even larger age spreads of tens of Myr, like in 
the Pleiades star cluster (Bclikov et al. 19981). The broad 
range of age spreads might indicate the existence of dif- 
ferent star formation scenarios for different star-forming 
environments. At the same time, discrepant results on 
the age spread in individual regions indicate that both 
observational and theoretical (e.g., Navlor 2009) biases 
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and the overall methodology used to assign ages to indi- 
vidual stars, might also be of importance. Observational 
difficulties include small number statistics, contamina- 
tion of samples by field stars, variable extinction and 
intrinsic infrared excess, unresolved binaries and insuffi- 
cient characterization of photometric uncertainties orig- 
inating in varying degrees of crowding. 

For our analysis, we have selected Westerlund 1 (Wd 1) 
and NGC 3603 YC, two of the most populous and mas- 
sive Galactic starb urst clusters (e.g.. iClark et al.l 120051 : 
iMelena et al.l [20081) . The clusters have half-mass radii 
of ~1 and 0.5 pc, respectively. They are composed of 
more than 10,000 stars each, a large fraction of which 
can be resolved individually via ground-based adaptive 
optics and Hubble Space Telescope (HST) observations 
from space. It is still unclear how such compact clus- 
ters have been formed and on wh ich time-scales. For 
NGC 3603 YC, IStohe et al.l (l200l suggest a single age 
of 1 Myr from isochrone fitting to the pre-main sequence 
(PMS) transition regio n and a sin gle burst of star for- 
mation, while iBcccari e t al.l (l2010( ) report a 10 Myr age 
spread for the PMS population, and therefore two dis- 
tinct episodes of star formation. F or Wd 1 recent studies 
give an age m the range 3-6 Myr (iBrandner et al.l l2008l: 
iGennaro et al.l[201lL iNeguerue la et alj 2010D. and^ n age 
spread of less than 1 Myr (^fegueruela etaD|2Qi3). The 
latter is based on spectral classification of Wdl's OB 
supergiant population. 

2. METHOD 

The first step of our method is a proper-motion se- 
lection of cluster members on the basis of mult i-epoch 
astrometric observations (e.g.. [Bedin et alll200lD . This 
enabled us to reject the majority of the contaminating 
field stars. Next, we apply Bayesian analysis to the 
photometry of cluster members with respect to theoret- 
ical isochrones to determine the probability distribution 
for the age of e ach member star, gi ven its photometric 
properties (e.g., IDa Rio et al.''2010'). We modified the 
Bayesian method of .Jorgenscn fc Lindcgren (,2005.) both 
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by taking into account the cluster membership probabil- 
ity and by adjusting the mass function (MF) by a com- 
pleteness factor. The posterior probability of the i-th 
star with magnitudes Ji^Ks^ to belong to the isochrone 
of age t is 



By multiplying the individual age distributions, we ob- 
tain the global probability function for the cluster's age t: 



L{t) = l[p{t\J,,Ks,). 



(9) 



p{t\J,,Ks,)^ p{J,,Ks^M,t)-aM\a,t)dM, (1) 



box 

where M is the initial stellar mass and £,{M\a,t) is 
the stellar MF with a slope a. The integration region 
is a "box", i.e., a rectangular area in color- magnitude 
(CM) space: Ks™„ < Ks < Ks^,,,iJ ~ Ks)rmn < 
(J — i^s) < ('^ — Ks)max- The limits of this area de- 
fine a mass range which is considered during integration. 
This box has been chosen in order to exclude field star 
contaminants based on their colors. To account for resid- 
ual field star contamination, the first multiplier inside the 
integral is defined as: 

p{J„Ks^M,t) = Pgood ■ p{Jr,Ks^M,t, cluster) 

+ {^- Pgood)-p{-h,Ks,\field), (2) 

where Pgood is the probability of a star to be a cluster 
member and (1 — Pgood) is the probability to be a field 
star. 

From the normalization conditions 



J p{Jt,Ksi\M,t,cluster)dJdKs = 1, 
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J p{J^,Ks,\field)dJdKs ^ I 
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we derive the probability that a cluster member or field 
star is found in a given position on the color-magnitude 
diagram (CMD) as 



p{Ji, Ksi\M, t, cluster) 



1 
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p{Ji,Ksi\ field) 
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(5) 
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A(A's)A(J-iirs)' 

where cr j. , axs ^-re the photometric uncertainties of 
star I and J{M,t), Ks{M,t) are the magnitudes of the 
theoretical isochrone, A{Ks) — Kg ~ Kg , A(J — 

Ks) = (J - Ks)max - (J - Ks)m^n■ Thc MF (iM\a,t) 

was modified by compl{M\t) to include source incom- 
pleteness and has the following form: 

^(M|a,i) = B ■ ■ compl{M\t), (7) 

where B was derived from the normalization condition 



^{M\a,t)dM = 1 



(8) 



and compl{M\t) from completeness simulations (see Sec- 
tion [S;!]). 



3. OBSERVATIONS AND DATA REDUCTION 
3.1. Observations of Westerlundl and NGC3603 YC 

Near-infrared adaptive optics observations of the cen- 
tral region of Wd 1 were carried out in 2003 April, using 
NACO at the Very Large Telescope (VLT). Ks observa- 
tions with a plate scale of 27 mas/pixel and a field of view 
(FOV) of 27" X 27" (corresponding to 0.5pcx0.5pc at 
4.0 kpc distance; (jGennaro et al.H20li[) ) were centered on 
R.A.(2000) = 16^47™06^5, decl.(2000) = -45°51'00". 
Ks frames with integration times of Iminute were co- 
added, resulting in a total integration time of 5 minutes. 

In 2010 August (epoch difference 7.3 yr), Wd 1 was ob- 
served with the HST Wide Field Camera 3 (WFC3/IR) 
in the F125W band with a plate scale of 130 mas/pixel. 
The final image consists of seven individual exposures 
with small (<10") offsets to compensate for bad pixels, 
with a total integration time of 2444 s. The overlapping 
FOV with NACO is k. 18" x 24". For a detailed descrip- 
tion of the full data set and data reduction, we refer to 
the work of M.Andersen et al. (2012, in preparation). 

Positions and ma gnitudes of th e stars were determined 
using DAOPHOT (jStetsonl 119871 ). Instrumental magni- 
tudes were calibrated against Two Micron All Sky Sur- 
vey, using suitable stars identifi ed on the NTT/SOFI 
data from the work by iGennaro et al.l (|201lD as sec- 
ondary photometric standards. To align the two data 
sets, identical bright stars with a high probability of 
being cluster members were identified. The geometric 
transformation of the NACO pixel coordinates to the 
WFC3 /IR system was based on fitting third-order poly- 
nomials with cross-terms in x and y. For bright stars 
with rriKs ^ 16-0 mag this resulted in an rms of RiO.036 
WFC3/IR pixel (Ri4.6mas) for the astrometric offset be- 
tween the WFC3/IR coordinates and the transformed 
NACO coordinates. 

Images of NGC 3603 YC were taken with HSTs Wide 
Field Planetary Camera 2 (WFPC2) with an image scale 
of 45.5 mas/pixel. The first epoch observations in the fil- 
ters F547M and F814W were separated by 10.15 yr from 
the second epoch observations in F555W and F814W fil- 
ters. Thc analysis was carried out for the core (<0.5pc) 
of NGC 3603 YC. The details of th e data reduction pro - 
cess have already been described bv lRochau et al.l (|201[K) . 

3.2. Completeness Correction 

As the matched HST /YLT data set for Wd 1 is lim- 
ited by the detection sensitivity in the J band (F125W), 
the completeness simulation was carried out only for the 
WFC3/IR data. Completeness simulations were done 
for the magnitude range mj from 14.0 to 23.0 mag. For 
each simulation, 25 artificial stars were added at random 
positions in the image, the image was then analyzed us- 
ing DAOPHOT, and the position and magnitude of the 
recovered stars recorded. This procedure was repeated 
1600 times, i.e. encompassing 40,000 artificial stars in to- 
tal. The recovery fractions were determined following the 
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steps outlined by iGennaro et al.l (|2011[ ). We calculated 
the average completeness for the WFC3/IR area overlap- 
ping with the NACO frame as a function of mj and fitted 
this by a Fermi function compl(M\t) = rnjjM^tf^ A^ — • 

e ^'2 +1 

Completeness simulations for NGC 3603 YC were done 
for F555W. Ten artificial stars were added at each of 50 
runs, in total 500 stars for each mp^^^w magnitude bin 
between 16.0 and 23.0 mag. We calculated the average 
completeness for _ff5T/WFPC2 image as a function of 
and fitted this by a Fermi function. 

3.3. Proper-motion Selection 

The main contaminants apparent in CMDs of the star- 
burst clusters are dwarf stars in the foreground and gi- 
ants in the background. Due to galactic rotation, their 
proper motion is different from that of cluster stars. Ve- 
locity disp ersions of disk and h alo stars are w50km/s to 
150km/s ()Navarro et al.|[20Tl . 

For Wd 1 (see Figure [1]) , our selection criterion for 
the astrometric residual of 5.8 mas corresponds to a 
proper motion of 0.8 mas/yr in the cluster rest frame 
(or 15.2 km/s at a distance of 4.0 kpc), which is higher 
than Wdl's interna l velocity dispersion of 2.1 km/s 
(jCottaar et al.l [20121 ). Such a selection provides an ef- 
fective discriminant between cluster members and field 
stars (see Figure [2^). 

For NGC 3603 YC, we used the result of IRochau et"all 
([2010), which is based on proper motions over an 
epoch difference of 10.15 years. The authors calcu- 
lated cluster members hip probabilities as described by 
iJones fc Walked ()1988[ ). and considered the stars with 
probabilities >90% as cluster members. 

4. CMD FOR WD 1 AND NGC 3603 YC 

The CMD for Wd 1 is presented in Figure [5^. For 
our further analysis we consider only the region with 
12.5 < mA's < 17.0 mag and 1.2 < mj — ttik^ < 2.9 mag 
(red box in Figure [2]), which comprises 41 stars with 
masses in the range from 0.5 to 11.5 Mq. Brighter 
stars were excluded because of saturation, fainter stars 
because of lower signal-to-noise ratio and hence larger 
photometric and astrometric uncertainties. The sam- 
ple includes main sequence (MS), Pre-MS, and tran- 
sition region stars, which have terminated their fuUy- 
convective Hayashi phase and are rapidly moving to- 
wards the MS. The CMD for NGC 3603 YC (Figure Eb) 
is derived from second epoch observations in F555W 
and F814W. For the age spread determination we se- 
lected the region with 16.5 < to_f555w < 21.5 mag and 
1.4 < mp555w — iTLpsiiW < 3.3 mag (red box), which 
comprises 228 stars with masses from 0.8 to 6.5 Mq. 
Overplotted is the best fitting isochrone assuming a par- 
ticular distance, for a solar metallicity Z=0.015, calcu- 
lated fro m the latest version of FRANEC evolutionary 
modelf0 ( Tognelh et alll2011[ ). adopting a mixing length 
value of ML= 1.68. We s upplement the FRA NEC mod- 
els with Padova models (jMarigo eTaDHOOl) for masses 
M > 7Mq. The FRANEC models have been trans- 
formed into the observational plane using spectra from 

^ FRANEC models in a wide range of masses and ages, for sev- 
eral chemical compositions and mixing length parameter, are avail- 
able at http://astro.df.unlpi.lt/stellar-models/ 



ATLAS9 model atmospheres (iCasteUi fc Kurucd |2004) . 
For the analysis, we used isochrones with 0.1 Myr spac- 
ing, covering an age range from 0.5 to 6 Myr. 

5. AGE LIKELIHOOD FOR WD 1 AND NGC 3603 YC 

The application of our method (Section ^ to Wd 1 
and NGC 3603 YC reveals a slight degeneracy between 
the cluster's distance and age. Since all cluster members 
should be at virtually the same distance, we can select 
a set of distances, and then analyze the resulting cluster 
age and age spread for each particular distance. 

In order to evaluate Pgood in Equation® for Wd 1 we 
estimated the density of stars in the CMD regions adja- 
cent to the red box (Figure [2l mj — rriK^ > 2.9 mag 
and mj — rriKs ^ 1-2 mag), where the stars are ap- 
parently non-cl uster members. The e xtinction value of 
Aks = 1.1 mag ([Brandner et al.ll2008D was assumed the 
same for all cluster members . The extinction law was 
taken from Ri eke fc Lebofskvl (|1985D . For the MF slope 
in Equation® we assumed a = 1.42 (where Salpeter 
slope would correspond to 2.3), as derived from near- 
infrared adaptive optics observations of the central region 
of Wd 1 (N.Kudryavtseva et al., 2012, in preparation). In 
order to quantify the photometric error a in Equation® , 
we used the results of the artificial star experiments to 
compare the known input magnitudes with output mag- 
nitudes recovered by DAOPHOT. A mo re detailed ex- 
lanat ion of the procedure is described in IGennaro et al.l 
20T1I) . The maximum photometric errors we got are 
0.05 mag in Kg and 0.18 mag in J for Wdl, and 0.17 mag 
in both bands of NGC 3603 YC. 

The L{t) function we derived from Equation® for 
Wdl at a distance modulus (DM) of 13.0 mag (4.0 kpc) 
is presented in Figure [3^. The full width at half maxi- 
mum (FWHM) of a Gaussian fitted to L{t) (red line in 
FigureOi) is 0.4 Myr. For 12.8 < DM < 13.2mag we got 
similar result. 

The L(t) function for NGC 3603 YC at DM=14.1mag 
is presented in Figure |4l As contamination by field 
stars was already significantly reduced by proper motion 
selection, Pgood in Equation® was estimated, assuming 
that cluster members are the stars with cluster member- 
ship probabilities >98%. The ext inction Ay = 4.9 mag is 
in a greement with t he res ults by iSung fc Besselll (|2004[ ) 
and IRochau et aTl (I2O1O0 a n d the relative extinction 
relations from 'Schlcgcl ct al." (1998j). We used the same 
a = 1.9 as in Stoltc ct al. (2006) for the MF slope. The 
maximum FWHM we derive for NGC 3603 YC at 13.9 < 
DM < 14.3mag is 0.1 Myr. 

6. BROADENING OF THE AGE LIKELIHOOD FUNCTION 

In the ideal case of a coeval population, which lies 
along an isochrone, and has no photometric errors, 
L{t) from Equation® would be a Dirac delta function. 
A number of observational and physical effects are 
potentially responsible for the L{t) broadening. In this 
section, we model these effects in order to estimate the 
true age spread. 

6.1. Photometric Error 

In order to quantify the broadening of L{t) due to 
solely photometric uncertainties we first generated a 
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number of cluster stars along an isochrone of a certain 
age in the CM space, and added random photometric 
errors. Random field stars were added to the data set 
with the same density as derived for real data. For this 
simulated data set, we applied the likelihood technique 
as described in Section [2] and got an artificial L{t) func- 
tion (see Figure Only due to photometric errors our 
simulation on 5.0 Myr population gave an L(t) broaden- 
ing in terms of FWHM equal to 0.25 Myr, hence the true 
age spread for Wdl should be even less than 0.4 Myr. 

6.2. Unresolved Binarity 

There is considerable observational evidence that bi- 
nary stars might consti tute a significant fraction of a 
cluster population (e.g.. iSharma et al]|2008 ). As an un- 
resolved binary combines the light of two stars, a binary 
will result in an offset in brightness by up to —0.75 mag in 
the CMD compared to a single star. For non-equal mass 
systems, there might also be an offset in color. Both cases 
lead to a broadening of the observed cluster sequence on 
the actual CMD. 

To test how the shape of L{t) is affected by unresolved 
binarity, we generated stars in the CM space in a sim- 
ilar way as described above for photometric errors, but 
with a 0.75 mag shift along the ordinate for 50% of the 
artificial points. We choose 50% as the mean value for 
the binarity fraction, as recent observations at least for 
massive population of Wd 1 revealed a high r ate (more 
than 40%) of binary stars ()Ritchie et al.ll2009D . We ap- 
plied the same Bayesian analysis to the simulated stars 
as to the real data to determine L{t). 

The normalized likelihood L{t) we obtained after com- 
bining the results of 30 simulations on binarity for a 
5.0 Myr population at DM 13.0 mag is shown in Fig- 
ure [3]:. It is clearly seen that binarity affects the shape 
of L{t) by adding a pronounced wing to the left from the 
main peak. Hence, a small shoulder toward younger ages 
from the L{t) maximum for Wdl (Figure [3^) could be 
caused by unresolved binarity. 

6.3. Ongoing Accretion 

For young stellar populations with ages < 10 Myr, on- 
going accretion and the accretion histor y of a star is also 
of imp ortance. As discussed in, e.g., iHosokawa et al.l 
Eom . low -mass objects, which gain mass through on- 
going accretion and a non-accreting PMS star of the 
same mass, arrive at different positions in a Hertzsprung- 
Russell diagram. 

For starburst clusters like Wd 1 and NGC 3603 YC, this 
effect seems to be strongly attenuated due to the pres- 
ence of very luminous and massive 0-type stars in the 
clusters. The fast winds and the ionizing radiation from 
these stars evaporate and remove circumstellar material 
around the low-mass cluster members, and clean out any 
remnant molecular gas in the cluster environment on 
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short timescales of a f ew lO^yr (e.g., [Adams "etall[200l 
iJohnstone et al.lll998D . This results in very little differ- 
ential extinction across the starburst cluster, enabling us 
to co nstrain stellar pro perties using broadband photom- 
etry (jStolte et al.ll2004[ ). 

6.4. Sensitivity to Age Spread 

In order to test the ability of our method to detect the 
real age spreads, we simulated a mixed cluster popula- 
tion with two different ages and assuming random pho- 
tometric errors. We repeated this procedure 30 times 
and calculated the averaged normalized likelihood L{t). 
The result for a mixed population with 70% stars of age 
5.0 Myr and 30% stars of age 4.5 Myr is shown in Fig- 
ure [3ji. Hence, a bimodal population with an age dif- 
ference >0.5 Myr would manifest itself in a prominent 
secondary peak, but the latter is not revealed in case of 
Wdl (see Figure |3i). 

7. SUMMARY 

Our analysis highlights the importance of a good mem- 
bership selection, background rejection and characteri- 
zation when analyzing crowded field data. We empha- 
size that the photometric component of our analysis only 
works in young, massive clusters and older clusters with 
little differential extinction and absence of ongoing accre- 
tion. Other environments require detailed spectroscopic 
analyses to establish precise effective temperatures and 
luminosities for individ ual stars, as has been exempli- 
fied in the case o f ONC ()Hillenbrandl[T997l ) and W3 Main 
(IBik et al.ll20T2h . 

Isochrones based on FRANEC and Padova evolution- 
ary models, and ATLAS9 atmospheric models provide 
a good match to observed cluster sequences in the age 
range from 1 to 5 Myr and mass range from 0.6 to UMq 
for solar metallicitics. Extending the sample of cluster 
members to lower masses could help to benchmark evo- 
lutionary tracks and atmospheric models for young, low- 
mass stars and quite possibly even brown dwarfs. 

Our analysis of the photometric sequence of cluster 
members yields an age spread of less than 0.4 Myr for 
the w5.0Myr old Wdl and less than 0.1 Myr for the 
w2.0Myr old NGC 3603 YC. This is a strong indica- 
tion that in both cases the clusters formed in a single 
event once a sufficient gas mass had been aggregated 
and compressed to overcome internal thermal, turbu- 
lent or magnetic support, and to initiate an avalanche- 
like star formation event. This finding seems to be 
in agreement with theoretical predictions that clusters 
with masses between 10^ and IO^Mq loose their resid- 
ual gas on a timescale s horter than their crossing times 
(|Baumgardt et al.ll2008h . For Wdl and NGC 3603 YC 
crossing times have been estim ated to be of the order of 
0.3 M yr (jBrandner et al.ll2008[ ) and 0.03 Myr (|Pang et al.l 
I2010D . respectively. 
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Fig. 1. — Proper motion diagram for Wdl. Stars with proper motions <0.8 mas/year are marked by green asterisks. 




Fig. 2. — Color- magnitude diagrams of Wd 1 (a) and NGC3603YC (b). Error bars indicate typical errors in color and magnitude. Red 
boxes show the regions that were taken for the analysis of the age spread. Proper-motion-selected stars for Wd 1 are marked by green 
asterisks. The CMD of NGC3603 YC includes only stars with cluster membership probabilities >90%. FRANEC-Padova isochrones (blue) 
are overplotted. For Wdl this is a 5.0 Myr isochrone at DM=13.0mag and for NGC3603 YC a 2.0 Myr isochrone at DM=14.1mag. 
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Fig. 3. — (a) Normalized L{t) for Wdl at DM=13.0 mag. The most probable age is 5.0 Myr. The red curve is a fitted Gaussian, (b) 
Normalized L(t) for the simulated 5.0 Myr population at DM=13.0mag. (c) Normalized L{t) for the simulated 5.0 Myr population at 
DM=13.0mag with a binarity fraction of 50%. (d) Normalized L{t) for the simulated population with 70% of the stars being 5.0 Myr old 
and 30% being 4.5 Myr old. 
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Fig. 4. — Normalized L{t) for NGC3603YC at DM=14.1mag. The most probable age is 2.0 Myr. The red curve is a fitted Gaussian 
function. 



